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We present a detailed study of the charmonium spectrum using anisotropic lattice QCD. We 
first derive a tree-level improved clover quark action on the anisotropic lattice for arbitrary quark 
mass by matching the Hamiltonian on the lattice and in the continuum. The heavy quark mass 
dependence of the improvement coefficients, i.e. the ratio of the hopping parameters ^ = KtfKs 
and the clover coefficients Ca,t, are examined at the tree level, and effects of the choice of the 
spatial Wilson parameter are discussed. We then compute the charmonium spectrum in the 
quenched approximation employing ^ = as/at = 3 anisotropic lattices. Simulations are made with 
the standard anisotropic gauge action and the anisotropic clover quark action with rs = 1 at four 
lattice spacings in the range as=0.07-0.2 fm. The clover coefficients Ca,t are estimated from tree-level 
tadpole improvement. On the other hand, for the ratio of the hopping parameters Q we adopt both 
the tree-level tadpole-improved value and a non-perturbative one. The latter employs the condition 
that the speed of light calculated from the meson energy-momentum relation be unity. We calculate 
the spectrum of S- and P-states and their excitations using both the pole and kinetic masses. 

We find that the combination of the pole mass and the tadpole-improved value of C, to yield the 
smoothest approach to the continuum limit, which we then adopt for the continuum extrapolation 
of the spectrum. The results largely depend on the scale input even in the continuum limit, showing 
a quenching effect. When the lattice spacing is determined from the IP — IS splitting, the deviation 
from the experimental value is estimated to be ~30% for the S-state hyperfine splitting and ~20% 
for the P-state fine structure. Our results are consistent with previous results at ^ = 2 obtained by 
Chen when the lattice spacing is determined from the Sommer scale rg. 

We also address the problem with the hyperhne splitting that different choices of the clover 
coefficients lead to disagreeing results in the continuum limit. Making a leading order analysis 
based on potential models we show that a large hyperfine splitting ~95 MeV obtained by Klassen 
with a different choice of the clover coefficients is likely an overestimate. 


I. INTRODUCTION 

Lattice study of heavy quark physics is indispensable for determining the standard model parameters such as the 
quark masses and CKM matrix elements, and for finding signals of new physics beyond it. Obtaining accurate results 
for heavy quark observables, however, is a non-trivial task. Since lattice spacings of order a « (2GeV)“^ currently 
accessible is comparable or even larger than the Compton wavelength of heavy quark given by 1 /rrig for charm and 
bottom, a naive lattice calculation with conventional fermion actions suffers from large uncontrolled systematic errors. 
For this reason, effective theory approaches for heavy quark have been pursued. 

One of the approaches is the lattice version of the non-relativistic QCD (NRQCD), which is applicable for a > l/rrig 
Since the expansion parameter of NRQCD is the quark velocity squared v^, lattice NRQCD works well for 
sufficiently heavy quark such as the the bottom (v^ ~ O.I), and the bottomonium spectrum |^-|^ and the bbg hybrid 
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spectrum have been studied successfully using lattice NRQCD. An serious constraint with the approach, 

however, is that the continuum limit cannot be taken due to the condition a > Ijrriq. Thus scaling violation from the 
gauge and light quark sectors should be sufficiently small. In practice it is often difficult to quantify the magnitude of 
systematic errors arising from this origin. Another difficulty is that there are a number of parameters in the NRQCD 
action which have to be determined. Since in the present calculations the tuning of parameters is made at the tree 
level (or tadpole improved tree level) of perturbation theory, the accuracy achieved is rather limited. 

Another approach for heavy quark uses a space-time asymmetric quark action, aiming to implement the 0{a) 
improvement for arbitrary quark mass [pd| . With appropriate parameter tunings, this action is unitarily equivalent 
to the NRQCD action im to higher order corrections for a > Ijmq, and goes over into the light quark Sheikholeslami- 
Wohlert (SW) action j^] for arriq <C 1. This approach has been originally proposed by the Fermilab group and the 
action is hence called the “Fermilab action”, whose first application is found in . Since the necessary tuning of mass- 
dependent parameters is in general difficult, in practice one uses the usual SW quark action even for a > Ijiriq, where 
the SW action is unitarily equivalent to NRQCD. This simplified approach, called the “non-relativistic interpretation” 
for the SW quark, has been widely used in current lattice simulations of heavy quark, such as the calculation of the 
B meson decay constant |TJ-|^. Toward the continuum limit a —> 0 the lattice action approaches the usual 0{a)- 
improved action and the systematic error becomes smaller as (aniq)^. However, the amq dependence at amg>l is 
quite non-linear, and it is not trivial how the systematic error could be controlled. 

Recently use of the anisotropic lattice for heavy quark simulations has been proposed |^,|^ as a possible alternative 
to solve the difficulties of the effective approach. On an anisotropic lattice, where the temporal lattice spacing at is 
smaller than the spatial one as, one can achieve atmq 1 while keeping agraq ^ 1. Therefore, using anisotropic 
lattices, one can reduce 0{{atmq)'^) (n = 1, 2, • • •) discretization errors while the computer cost is much less than that 
needed for the isotropic lattice at the same at- Naively it is expected that the reduction of 0{{atmq)'^) errors entails 
the reduction of most of discretization errors due to large quark mass, since the on-shell condition ensures that the 
large energy scale flows only into the temporal direction as far as one considers the static particle, with zero or small 
spatial momentum. If such a naive expectation is correct, the discretization error is controlled by a small parameter 
atmq as it is for light quarks, and one can achieve even better accuracy by taking a continuum limit. However, it is 
not obvious that one can eliminate all 0((asTOq)"') errors at the quantum level, even if it is possible at the tree level. 

Another advantage of the anisotropic lattice, which is more practical, is that a finer temporal resolution allows us 
to determine large masses more accurately. This has been already demonstrated in simulations of the glueball [pO 21 
and the hybrid meson . 

Klassen calculated the charmonium spectrum in the quenched approximation, employing lattices with the ratio of 
the temporal and spatial lattice spacings ^ = as/at = 2 and 3, as a feasibility study of the anisotropic approach |]l8| , |l^ . 
He tuned the ratio of the temporal and spatial hopping parameters (/ = Kt/Kg non-perturbatively by demanding the 
relativistic dispersion relation for mesons. For the spatial clover coefficient Cg, he adopted two choices, the tree level 
tadpole improved value correct for any mass {atmq > 0) and that correct only in the massless {atmq = 0) limit, in 
order to make a comparison. He mainly studied the spin splitting of the spectrum, and obtained an unexpected result 
that two different choices of the clover coefficients lead to two different values of the S-state hyperfine splitting even in 
the continuum limit |r^ , [T^ . The continuum limit is of course unique, and clearly, at least one of the two continuum 
extrapolations is misleading. Since the hyperfine splitting is sensitive to the clover coefficients, it is plausible that the 
disagreement is due to a large discretization error arising from the choice of the clover coefficients. In an unpublished 
paper |^, he pointed out the possibility that the 0{{^atmq)'^) = 0{{asmq)'^) errors still remain with his choice of 
the parameters, which we review in the next section. A similar statement can be found in some recent studies 
In fact, he adopted rather coarse lattice spacings Og ~ 0.17—0.30 fm where Ogm^ ^ 1. It is then questionable whether 
the reliable continuum extrapolation is performed at such coarse lattice spacings. 

Using the same anisotropic approach as Klassen, Chen has recently calculated the quenched charmonium spectrum 
(Q. She employed ^ = 2 and finer (og ~ 0.10—0.25 fm) lattices, and adopted the tree level tadpole improved 
clover coefficient Cg correct for any mass, which is expected to be better than the other choice that is correct only 
in the massless limit. She computed not only the ground state masses but also the first excited state masses, and 
extrapolated them to the continuum limit. Her results at ^ = 2 are consistent with Klassen’s results at ^ = 2 and 3 
with the same choice of the clover coefficients. 

Since Chen’s calculation was performed only at ^ = 2, similar calculations at different values of ^ using fine lattices 
are needed to check the reliability of the continuum limit from the anisotropic approach. In addition, the complete 
P-state fine structure has not yet obtained in this approach so far, since the mass of ^P 2 {Xc 2 ) state has not been 
measured in previous studies. 

In this work, we present a detailed study of the charmonium spectrum from the anisotropic lattice QCD. We perform 
simulations in the quenched approximation at ^ = 3, employing fine lattice spacings in the range Og = 0.07—0.2 fm. 
We attempt to determine the ground state masses of all the S- and P-states (including ^^ 2 ) as well as their first 
excited state masses. To estimate the systematic errors accurately, we adopt both the tree level tadpole improved 
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value and non-perturbative one for C, and both the pole mass and kinetic mass for Miat(l<§) which is tuned to the 
experimental value. We focus on the lattice spacing dependence and continuum limit of the mass splittings. We 
compare our results with the previous anisotropic results by Klassen and Chen to check the consistency, and with 
experimental values p5|| to estimate the quenching effect. 

In addition, to understand the discrepancy of the hyperfine splitting mentioned above, we make a leading order 
analysis using the potential model. To examine the effect of clover coefficients, we estimate the hyperfine splitting at 
leading order. Comparing the leading order estimates with numerical results for the hyperfine splitting, we attempt 
to find a probable solution for this problem. Our preliminary results are already reported in Refs. |26|j27| . 

This paper is organized as follows. In Section we summarize and discuss the theoretical aspect of the anisotropic 
lattice QCD. In Section III, we give details of our simulation. Our results for the charmonium spectra are shown in 
Section IV, where we attempt to take the continuum limit and estimate the quenching effect. We address the problem 
of the discrepancy of the hyperfine splitting and study the effect of clover coefficients in Section Section VI is 
devoted to our conclusions. 


II. ANISOTROPIC LATTICE QCD ACTION 

In this section we first define the anisotropic lattice action used in this work and fix notations. We then derive the 
tree level values of bare parameters in our massive quark action, and discuss effects of the anisotropy. Although it was 
already discussed in earlier papers we briefly describe the outline of derivations in order to be self-contained. 

We also consider the tadpole improvement of bare parameters and see how tree level values are modified. 


A. Anisotropic gauge action 


In this work, we use the standard Wilson gauge action defined on an anisotropic lattice: 




1 ^ 

^ x,s>s' x,s 


Pst{x)) 


( 1 ) 


where f3 = 6 / 5 ^ is the gauge coupling, and Pss'{x) and Pst{x) are the spatial and temporal plaquettes with Pfj,u{x) = 
^Re Tr U^y{x). The anisotropy is introduced by the parameter ^0 and we call this the “bare anisotropy”. We 
denote spatial and temporal lattice spacings as Og and at and define the “renormalized anisotropy” ^ = asjat- We 
have ^ = ^0 at the tree level, and the ^ = C(Co, /3) at finite (3 can be determined non-perturbatively by Wilson loop 
matching [p 8 |-^. In numerical simulations, there are two methods for anisotropy tuning, either varying ^0 to keep 
^ constant or vice versa. Since the former is more convenient for keeping the physical size constant, and easier for 
performing the continuum extrapolation, we adopt it in this work. 


B. Anisotropic quark action 


For the quark action, we employ the space-time asymmetric clover quark action on an anisotropic lattice proposed 
in Ref. jllQ: 


X 


( 2 ) 


Q = Too -f r'oWoVo 


Co 




woy^croi.Foi(a;) + 7^ V' 

x,i<j 


>/. 

x.i 


( 3 ) 


where vq = 1 and mg = atVfiqQ is the bare quark mass, and and F^i, = a^a^F^i, with (oo,ai) = 

{at, Us). The Wilson operator IF^ is defined by 



(m = 0,I,2,3) 


with the Wilson coefficients {ro,ri) = {rt,rs) and 
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( 4 ) 









( 5 ) 


- 2a, 
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— 2 ^fj.,x'4’x+jl + x-jl4’x — jl 2 ,'iIJx 


( 6 ) 


For the field tensor we adopt the standard clover leaf definition. Note that, in Eq. (^), the factors in front of spatial 
Wilson and clover operators include rather than This is merely a convention and there is no deep theoretical 
reason. This action is essentially the same as the one employed by Klassen and Chen |^. In Chen’s work, 
however, vq was a tuning parameter with z/ = 1 fixed. The two parameterizations are related to each other by a field 
rescaling ipx ^ '4’xl'Jv. Therefore {toq, z^o, w, wo}|3 corresponds to {mo/zz, 1/z/, w/z/, in our convention. Among 

these six parameters {mo, v, Ts, rt,u!, wq}, at least one is redundant, so that we take r* as a redundant parameter and 
use it to remove the fermion doublers. Although Vg may not be taken arbitrary in the 0{a) improved anisotropic 
quark action for the matrix elements, it can be taken arbitrary for the hadron mass calculation. Therefore we 
always set rt = 1 and leave free in this work. The remaining parameters {mo, z^, are used to tune the quark 

mass and reduce the lattice discretization error. 

For convenience in numerical simulations, we also present the quark action in a different form. Rescaling the fields 
il^x, the quark action can be transformed into a form given by 


S'f = '^{^x'^l^x - - 'yo)Uo,x'ipx+d + '^x{l+ 1 o)uI^^_q^x- 6\ 

X 

- KsY^[4>x{rs - 'yi)U^,x'4’x+i + 4>x{rs + li)Ul^_i^x-i\} 

i 

^ ^ '^x^ij ^ij ~F ^ ^ '^x'^0i^0i(,^')'4^x j 

x,i<j x,i 


(7) 


where Kg^t and are the spatial and temporal hopping parameters and the clover coefficients, respectively. The 
hopping parameters ATs.t are related to the bare quark mass mo = Utm^o through 


atm,o = l/{2Kt) - 3r«/C - 1, C = Kt/Kg. (8) 

The form Eq. (Q) on the anisotropic lattice is the same as that on the isotropic lattice in Ref. ||^ . Note however that 
Ref. uses the inverse of our definition for (. We refer to their definition as Cf = Kg/Kt = VC- Using Eq. (^) one 
can convert {mgo,C} to {Kg,Kt\. In our convention, the relation between {z^,a;,a;o} and {C, Cs,Ct} is given by 


C = Cg=ujlv, Ct = ioU}o/v 


(9) 


or equivalently. 


= ‘?o/C> ^ azo = ctvHo¬ 


rn 


Following Ref. 0, we call the quark action Eq. (^) as the “mass form” and Eq. (i as the “hopping parameter 
form”. 


C. Tree level tuning of bare parameters for arbitrary mass 


To derive the tree level value of bare parameters, we follow the Fermilab method and calculate th e la ttice Hamil¬ 
tonian After some al gebra (see appendix ^ for details), we obtain the lattice Hamiltonian Eq. (M). Using the 
FWT transformation Eq. (A17), we then transform it to the non-relativistic form, in which the upper components of 
the Dirac spinor completely decouple from the lower ones {i.e. eliminate y-D and a-E). The transformed Hamiltonian 
is given by 


—H'^ = 
at 


mi -b 7 oAo - 


D2 

2m2 


iS • B [7 • D, 7 • E] 

2ms 


■ 7o- 


8m I; 



( 11 ) 


^More precisely, Chen used the language {mo, z^t, Cgwi Usw} instead of {mo,uo,uJ,LOo} 
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with 


OtTOi 

1 

atm2 

1 

atniB 

1 

{atmEY 


= log(l + TOo), 


+ 


^sCf 


mo(2 + TOo) 1 + mo’ 


2Cf' 


+ 


CsCf 


mo (2 + TOo) 1 + mo’ 

(1 + moY 


= 4C^ 


mo (2 + m^Y 


+ (ct - 1 ) 


mo(2 + mo)_ 


( 12 ) 

(13) 

(14) 

(15) 


where Cp, T'g and c), are defined in Eq. (AS). The S • B term gives the leading order contribution to the hyperfine 
splitting, while [7 • D, 7 • E] term yields the fine structure splitting. 

The matching condition = H^r + 0{al) is equivalent to 


mi = m2 = ruB = niE = rriq- 

This yields the tree level value of bare parameters for massive quark: 


^oCf = v = 


/ ^or^mo(2 + mo) ^ 

V 4(1 +mo) ) 


mo(2 + mp) 
21 og(l + mo) 


Co^’s™o(2 + mp) 

4(1 +mo) 


(16) 


(17) 


{lU = Tsiy), 


(18) 


^ {^oCfY - 1 , ^o^sCf (^or,)^mo (2 + mo) 

* mo(2 + mo) 1 + mo 4(1 + mo)^ 

We note that Cg is independent of the quark mass, while v and ct have complicated mass dependences. The term 
^omo — asmqo seems to exist in Eq. © and ©)■ To see this explicitly, we expand e and ct in mo. This gives 

= 1 + ^(1 - iors)mo + ^(-1 + ^^ors + 3(^o»'s)^)rno + O(mo), ( 20 ) 


Ct = ^ ~ ^^orsY)mo + 0{ml). ( 21 ) 

The ttsmqo term, which is 0(1) for heavy quarks at currently accessible lattice spacings of « 2 GeV, appears in 
E and Ct even at the tree level. Since ^o™o = ctgmqQ is always multiplied with the spatial Wilson coefficient in 
Eqs. (p^) and ([^), one can eliminate the asiriqo term at the tree level by choosing 

Ts = 1 /^ 0 . ( 22 ) 

However, this choice has the disadvantage that the mass splitting between unphysical doubler states and the physical 
state decreases as ^0 increases. Moreover, the hopping terms in the quark action are no longer proportional to the 
1 + 7 ^ projection operators. It is also doubtful that, beyond the tree level, the Ogm^o term can be still eliminated by 
this choice. 

If one adopts the conventional choice 


rs = 1, (23) 

the Osmgo term remains, but the unphysical doubler states decouple. This choice also has the practical merit that 
the quark action has the full projection property, so that the coding is easier and the computational cost is lower. 

The tree-level full mass dependences of e and Ct for Vg = l/^o and rg = 1 are shown in Figs. ^ and |[ In order to 
compare at the same Ug, we choose miOg as the horizontal axis instead of mi at where mi is the pole mass. Since 
aj^>l GeV and mi < mbottom ^ 4.5 GeV in current typical simulations, we plot results for miOg < 4. 

For Tg = 1/Co shown in Fig. |^, both e and Ct are monotonic functions in mass, and they converge to their massless 
values as Co increases at any fixed values of miCg. Hence, the agiriqo dependence can be controlled by increasing Co- At 
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^ = 100 the mass dependences of v and ct completely disappear with the cost that the physical and unphysical states 
are almost degenerate. In actual simulations with = l/^Oi taking 2 < oo to decouple unphysical doublers, 

one is allowed to use the massless values for v and Ct, since their mass dependences are monotonic and very weak. In 
this case mass dependent parameter tuning can be avoided even at agrrio ^ 1. 

For Tg = I, on the other hand, the mass dependences of v and ct are complicated and non-negligible even for large 
^ 0 - Indeed i/ and c* do not converge to their massless values as increases at fixed miCs, as shown in Fig. The 
deviation from the massless values at = 2 is smaller than the one at = 1, but it becomes larger again as 
increases. Therefore, taking = 2-3 in simulations with = 1, one needs to perform a mass dependent parameter 
tuning. 

For both choices of r^, it is better to use a moderate value of rather than excessively large values. In our 

numerical study of the charmonium spectra, we adopt the choice = 1, and make a mass dependent parameter 
tuning, due to the practical reasons mentioned above. 

Finally we show the tree level value of the parameters in the massless limit. By taking atrriqQ —> 0 in Eqs. 0)-®, 
one obtains 


in the mass form, or 


j/ = 1, w = rs, LOO 


1 + ^prs 
2^0 


C = Co, Cs = rs, Ct 


1 + Co^s 
2 


(24) 


(25) 


in the hopping parameter form. Note that there is an ambiguity in the tree level value of as/at, since = C ^-t the 
tree level but fo C in the simulation. Fortunately, this ambiguity almost disappears after the tadpole improvement, 
as shown in the next subsection. 


D. Tadpole improvement 

In this section we apply the tadpole improvement ||^ to the parameters of the anisotropic lattice action at tree level 
in order to partially include higher order corrections. One first rewrites the lattice action using a more continuum-like 
link variable Ui^ = Uifi/ug^t, where Ug^t = {Uifi) is the expectation value of the spatial or temporal link variable, i.e. 
one replaces 

Uifi Us,tUifi, (26) 

and then repeats the tree-level calculations. We will show below how the tree-level values of bare parameters are 
modified. 


1. Gauge action 


By the replacement of Eq. (p^, the anisotropic gauge action Eq. (Rl) becomes 


-E 


6 


Pss' + Co Pst + const, independent of Up 


where = iRe Tr and and are given by 




Ulut y'{Pss') {Pst) 


Co = -Co ^ 

Us 


(Pst) 

{Pss') 


Co- 


(27) 


(28) 


Requiring space-time symmetry for the action Eq. ( 0) in the classical limit, one obtains the tree-level tadpole- 
improved value of the anisotropy (denoted by an index ‘TI’), 


C"*"^ = Co = {ut/us)Co- 


(29) 


In practice in Eq. (29) agrees with the renormalized anisotropy C within a few % accuracy at ~ 1. Therefore one 
can replace the factor {ut/us)Co by C in the following equations. This simplifies the tree level expression. Moreover, 
the arbitrariness for the choice of anisotropy disappears. 
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2. Fermion action 


When the fermion action is rewritten in terms of Ui and {/q instead of Ui and t/p, the action keeps the same form 
with 


Ks = UsKs, Kt = UtKt, 

(30) 

Cs=MgCs, Ct=UsU^Ct. 

(31) 

Then Q = Kt/Kg and the bare quark mass atrriqg = l/2Kt - (Id- 3rs/C) are modified to 


c = Kt/Ks = {ut/us)C, 

(32) 

atrhqg = - (1 + 3rs/C) 

ZKt 

Ut Ut Ut 

(33) 

Using parameters with the tilde, one can repeat the derivation in the previous subsection, 
obtains 

For a massless quark, one 


1 + 1 + 


C = Co - C. Cs = Ts, Ct = 

z z 

Therefore, tadpole-improved (TI) tree-level estimates are 

= {Us/ut)io = Co, 

which indicates that non-perturbative C at fhqo ~ 0 is closer to Co than to C, and 

1 1 + {ut/us)^ors 1 1 + C^s 


„TI 


- 


UsUl 


UsUf 2 


(34) 


(35) 


(36) 


As can be seen in Eqs. (^^ and (|^), the tadpole improvement eliminates the uncertainty of choice of anisotropy {i.e. 
whether to chose Co or ^) at tree level. Converting to the {z/, WjWo} convention, one obtains 


= 1 , uj'-^ = = 


.TI 


TI 


1 1 - 1 - {ut/us)^ors 

ulut 2{ut/Us)^0 


(37) 


Note that is normalized to 1 since i/ equals Co/C and not C/C, hence the former definition is practically more 
convenient than the latter one. Note also that tadpole factors in cj^ and are different because cog equals Ctv 
and not Ctv 

Similarly, for massive quarks, tadpole-improved tree-level estimates become 


1 /C^i = - 

Us 


J 1 /'rsrfig{2 + rhg)^ 

rhg{2 + rhg) 

rsrhg{2 -|- rhg) \ 

|VV 4(1-two) , 

1 2{ut/us)‘^^g\og{l + rhg) 

4(l-|-rfio) j 


(38) 


with z/Ti = ^o/C , and 


c? = 


(39) 


„TI _ 
— 


1 (z/T1)2 _ 1 

su} I mo(2 -I- Too) 


Us J 1 -I- Too 


(Cor^)^mo(2 -l-mo) 

Us) 4(1-|-too)^ 


(40) 


where rfig = atrhqg. 
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III. SIMULATIONS 


We proceed to calculate the charmonium spectrum in the quenched approximation as our first numerical study 
using the anisotropic lattice. In this section we describe the computational details of our quenched charmonium 
calculation. 


A. Choice of simulation parameters 


For the gauge sector, we use the anisotropic Wilson gauge action given in Eq. (||). Throughout this paper, we 
employ ^ = 3, where ^ is the renormalized anisotropy. In order to achieve ^ = 3, we tune the bare anisotropy 
using the parameterization oi rj = given by Klassen 


viPA) = 1 + 1-7 —r~ - 2 5 > 

V 4 / 6 1 + aog^ 


where qq = —0.77810, ai = —0.55055 and 


miO = 


1.002503^3 0.39100^2 1 . 47130 ^ _ 0.19231 

^3 + 0.26287^2 1.59008^ - 0.18224 ' 


(41) 


(42) 


We perform simulations in the quenched approximation, at four values of gauge coupling (3 = 5.70, 5,90, 6.10 and 
6.35. These couplings correspond to as = 0.07—0.2 fm and atmcharm = 0.16—0.48 for mcharm = 1-4 GeV. The spatial 
lattice size L is chosen so that the physical box size is about 1.6 fm, while the temporal lattice size T is always set to 
be T = 2^L = QL. 

For the charm quark, we use the anisotropic clover quark action Eq. ( 0 ) with the conventional choice of the spatial 


Wilson coefficient, = 1, as mentioned in Sec.IIC, We take two values for the bare quark mass uiq = (toJ,TOq) at 


each f3 in order to interpolate (or extrapolate) results in mo to the charm quark mass The charm quark mass 

^charm jg gxed from the experimental value of the spin averaged IS' meson mass. In this procedure, we use both the 
pole mass Mpoie and kinetic mass Mkin for the IS meson. Eor C, the ratio of the hopping parameters, we adopt both 
the tree-level tadpole-improved value ^^nd a non-perturbative value determined from the meson dispersion 
relation. We describe our method of tuning C in detail in Sec. HI C| . Eor the spatial clover coefficient Cg, we employ the 
tree-level tadpole-improved value for massive quarks Eq. Note that Cg has no mass dependence at the tree level. 
On the other hand, we adopt the tree-level tadpole-improved value in the massless limit Eq. (^ for the temporal 
clo ver co efficients c*. We discuss possible systematic errors arising from our choice of the parameters C and Cg^t in 
Sec.Ill E. The tadpole factors Ug^t in Eqs. (p^) and (|3^) are estimated by the mean plaquette prescription: 


= {Pgg')^/\ ut = 1 . 


(43) 


If we adopted the alternative definition ut = {{Pgt)/{Pss'))^^^ instead, Ut would be greater than 1. We use ^ instead 
of {ut/ug)^o in Eq. (||). 

Gauge configurations are generated by a 5-hit pseudo heat bath update supplemented by four over-relaxation 
steps. These configurations are then fixed to the Goulomb gauge at every 100-400 sweeps. On each gauge fixed 
configuration, we invert the quark matrix by the BiGGStab algorithm to obtain the quark propagator. We always 
perform the iteration of the BiGGStab inverter by T times, where T is the temporal lattice size. By changing the 
stopping condition for the quark propagator, we have checked that this criterion is sufficient to achieve the desired 
numerical accuracy. We accumulate 400-1000 configurations for hadronic measurements. 

Our simulation parameters are compiled in Tables | and ||. In Table III, we compare some of the parameters used 
in our simulation (labeled by “set A”) with those in the previous studies by Klassen (“set B” and “set D”) [pfr 19 
and by Chen (“set C”) Q for later references. 


B. Meson operators 

In this work, we calculate all of S- and P-state meson masses of charmonia, namely ^SqAc), ^Si{J/tj}), ^Pi{hc), 
^PoiXco), ^PiiXci) and ^P 2 {Xc 2 )- For this computation, we measure the correlation function of the operators which 
have the same quantum number as one of above particles. In Table 0 we give the operators for the S- and P-state 
mesons. There are two types of operators, those of the form ^F-f/; and of i/jrA^, where F represents a combination of 
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7 matrices and A the spatial lattice derivative. We call them the F operator and the FA operator, respectively. The 
latter appears only for the P-state mesons. Note that there are two lattice representations for the ^P 2 state (E and 
T representations) due to breaking of rotational symmetry. 

We measure the correlation functions of the F operators 

Qtlte(0 = • E ^-o.orV'yo.o/X-zo/xo-yo). (44) 

X yo.zo 

where /® is a source smearing function, and we always adopt a point sink. We employ the point source (s = 0) with 
= (5x,o and an exponentially smeared source (s = 1) with where As and Bg are smearing 

parameters. Therefore we have three source combinations, ss' = 00, 01 and 11, for the F operators. The smearing 
parameters Ag and Bg at each /3 are chosen so that the effective mass of the 15 meson for ss' = 01 has a wide plateau. 
To obtain the correlation functions of the FA operators, we measure 

(45) 

X yo,zo 

where Ai'i/jx,* = V'x+i t ~ V’x-i t ^4ie discretized derivative at the sink, and we employ a smeared derivative source 
(s = 2) given by 


- Age-^‘‘^^-^\ {i = 1, 2, 3) 


(46) 


with Ag and Bg the same as those for s = 1. For the ^Pq state, for example, we calculate Clp^ = ^ujj 

Fi = 7 i. For the FA operators, we have two source combinations, ss' = 02 and 12. In total, S-state mesons have 
ss' = 00, 01 and 11 source combinations, and P-state mesons have 00, 01, 11, 02 and 12 source combinations except 
for ^P 2 - Since there is no F operator for ^P 2 , it has only 02 and 12 source combinations. 

To calculate the dispersion relation of S-state mesons, we measure correlation functions for four lowest non-zero 
momenta. 


OsP = 


(27r/i)x {(1,0,0), (1,1,0), (1,1,1), (2,0,0)}, 


(47) 


in addition to those at rest. Correlation functions with the same value of |p| but different orientations are averaged 
to increase the statistics. 


C. Tuning bare quark mass mo and fermion anisotropy (( 


Let us describe our method of tuning ( and uiq in detail. We determine the input parameters mQ{= ml,mQ) and 
C(= C^^) as follows. First we fix C = ^ = 3 and choose and Wq where the IS meson mass roughly agrees with 

the experimental value. Then we determine both the tree-level tadpole-improved value and the nonperturbative 
value at itiq = wj and m§. 

To obtain at fixed mo , we use Eqs. ( p^ and (^^. We replace the factor ut/ug in Eq. ( |^ ) with using 

Eq. (^^. On the other hand, is obtained by demanding that the relativistic dispersion relation is restored at 
small momenta for the IS meson. The dispersion relation on a lattice is given by 


E{pf = E(Qf + c2p2 + 0(ay) 


(48) 

(49) 


where c is called the ‘speed of light’, and Mpoie and Mkin are the pole and kinetic masses of the IS meson. Throughout 
this paper, a capital letter M denotes the meson mass, while a small one m the quark mass. Generally c is not equal 
to one due to lattice artifacts. We extract the speed of light c by fitting E{pY linearly in p^ for three or four lowest 
momenta, since the linearity of E{p)^ in p^ is well satisfied. We identify with a point where c = 1 or equivalently 
Mpoie = Mkin for the IS meson. To determine we perform preparatory simulations and calculate c for ^ = 2.8, 
3.0 and 3.2 at mo = mj and mp using 100-200 gauge configurations. Then we find C = where c = 1, from an 
interpolation of C. As shown in Table |l|, the speed of light c at is indeed equal to 1 within 1%, which is roughly 
the size of the statistical error. 
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Production runs for the charmonium spectrum described in Sec .|lII A are performed at rriQ = (mQjTOp) and C = 
(^Ti^ ^NP) Accidentally, for /3 = 5.90 and 6.35, holds within our numerical accuracy, so we use 

the same data for the analysis at these (3. 

Finally we linearly interpolate or extrapolate results at toq = (mj, tuq) to those at mo = mQ**”™, with fixed C (= 
or As already mentioned, we identify mo*'“™ with a point where the spin-averaged IS meson mass Mis^ilS) in 

units of a physical quantity Qiat is equal to the corresponding experimental value: 


Mat (15) 

l^lat 


Mexp(15) 


Q 


exp 


(50) 


with Mexp(15) = 3067.6 MeV for charmonium. In this work, we adopt the Sommer scale tq and the spin-averaged 
mass splittings AM {IP — 15) = M (IP) — M(15) and AM (25 — 15) = M (25) — M(15) as the scale quantity Q. The 
spin-averaged masses are defined by 

M{nS) = [3M{n^Si) + M(nMo)]/4, (51) 

M{nP) = [3M(nMi) -h 5M{n^P2) + 3M{n^Pi) + M(nMo)]/12 (52) 

with n{= 1,2,- • •) the radial quantum number. The experimental values of the mass splittings AM {IP — 15) and 
AM(25 — 15) are 457.9 MeV and 595.4 MeV, respectively. The experimental values of tq is not known, and we use a 

phenomenological estimate rp = 0.50 fm. For the definition of the lattice meson mass Miat in Eq. (^0|), we have two 

choices in the case of (^ = one is the pole mass Mpoie and the other is the kinetic mass Mkin- On the other hand, 
in the case of C = Mpoie = Mki^ should hold by definition. In practice, there can be small deviations due to the 
statistical error. Therefore we have 4(=2x2) choices for (Miat,C) in total. 


D. Mass fitting 

From meson correlation functions we extract the meson mass (energy) by standard fitting with a multi-hyperbolic- 
cosine ansatz (termed rifit-cosh fit below) 

^fit 1 rji 

C'sTate(i) = A'"'cosh[(--t)Mi], (53) 

i=0 

where ss' represents the source combination ( 00 , 01 , etc.), t is the time separation from the source, and nut is the 
number of states included in the fit. 

We determine the mass of the ground state and the first radial excited state for each particle, and the mass splittings 
such as AM {IP — 15) and AM(25 — 15), from a 2-cosh fit using several correlation functions with different source 
combinations simultaneously. Here we use the correlation functions of ss' = 00, 01 and 11 sources for S-states, while 
00, 11, 02 and 12 sources are used for P-states except for ^P 2 - For ^P 2 , we use the correlation functions of 02 and 
12 sources. The 2-cosh fit for each S-state always gives the ground state mass consistent with that from the 1-cosh 
fit. On the other hand, for the P-state, the 2-cosh fit is preferred over the 1-cosh fit because the IP mass from the 
1 -cosh fit using the correlation function of 11 and 12 sources occasionally disagrees by a few cr, due to excited state 
contaminations. To determine the mass of the first excited state accurately, it is better to adopt results from the 
3-cosh fit. However, we do not perform the 3-cosh fit systematically because of the instability of it, and adopt results 
from the 2-cosh fit for the first excited state mass. This may cause an overestimation of the first excited state mass 
due to a contamination from higher excited states. 

To determine the spin-averaged 15 mass and the 15 energy at p 7 ^ 0, and the spin mass splittings such as AM(1^5i — 
1^5o) and AM(l^Pi — I^Pq), we perform a 1 -cosh fit (ngt = 1 ) using the source combination which gives the widest 
plateau in the effective mass. We use the 01 source for the S-state, and the 12 source for the P-state. We always 
check that the spin mass splitting from a simultaneous 2 -cosh fit mentioned above agrees with that from the 1 -cosh 
fit within 1-2 a. We also check that the splitting AM(l^Pi — l^Po) from a 1-cosh fit using the 11 source agrees with 
that using the 12 source. 

In these analyses, we perform both the uncorrelated fit, and the correlated fit which takes account of the correlation 
between different time slices and different sources. The uncorrelated fit is always stable and gives /Ndf<QA {Q ~ 
1). The correlated fit with 1-cosh ansatz is also stable and produces results consistent with those from the uncorrelated 
fit. However, the correlated 2-cosh fit is often unstable, either failing to invert the covariance matrix or giving large 

/Ndf ^ 1 even if it converges. Therefore we adopt the uncorrelated fit for our final analysis. 
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The fitting range [tmin, tmax] for the final analysis is determined as follows. From an inspection of the effective mass 
plot, we determine tmax which roughly has the same physical length independent of j3. We repeat the 1- and 2-cosh 
fits for each /3, varying tmin with fixed tmax, and find a range of tmin where the ground state mass and the first excited 
state mass (for 2-cosh fit) are stable against tmin- We also check that it has reasonable value of x^/^df- The final 
tmin is then chosen from the region accepted above so that its physical length is roughly equal independent of B. 

Typical examples of the effective mass plot and tmin-dependence of the htted mass are shown in Figs. and 
in Fig. H, respectively. Our final fitting ranges are summarized in Table Statistical errors of masses and mass 
splittings are estimated by the jack-knife method. The typical bin size dependences of jack-knife errors for the ground 
state masses are shown in Figs. ^ and ^ We always adopt a bin size of 10 configurations, i.e. 1000-4000 sweeps. 


E. Scaling violation and the continuum limit 

We discuss scaling violation for our action and how the results at finite Og are extrapolated to the continuum limit 
as 0. Since we use the anisotropic Wilson gauge action with nonperturbatively tuned the scaling violation from 
the gauge sector starts at 0((asAQCD)^)- 

For the quark sector, we use the anisotropic clover quark action with tadpole-improved clover coefficients Cg,*, 
and either the tadpole-improved value or nonperturbative value for Since we adopt the tree-level tadpole- 
improved value of Cg for massive {agmq > 0) quarks, the scaling violation arising from the choice of Cg is 0((asAQCD)^) 
and O(aasAQCD)- On the other hand, for ct, we adopt the tree-level tadpole-improved value correct only in the 
massless (ogTOg = 0) limit, which generates an additional 0(agAQCD • OgTOq) = 0(0^AgcD^T^g) error. Recall that the 
UsTTig (not only atrug) dependence of the parameter remains with our choice of the spatial Wilson coefficient rg = 1 at 
the tree level, as discussed in Sec.|i[ In the case of C = therefore, the scaling violations are 0((asAQCD)^) and 
0(agAqcD^Tig) at leading order, and O(aagAQCD) at next-to-leading order. The size of these errors are estimated to 
be 0((agAQCD)^) = 7—1%, O(asAQCDWg) = 37—4% and O(aOsAQCD) = 4—1% for (3 = 5.70—6.35 corresponding to 
« 1.0—2.8 GeV. Here we took Aqcd = 250 MeV (~ A^“°) and mq = 1.4 GeV (~ Wcharm), and the renormalized 
coupling constant a is estimated from Eq. (^). It is expected that the O(aagAQCD) errors are largely eliminated by 
the tadpole improvement. 

When the tree level tadpole improved value is used instead of we have additional 0{a) and 0{aasmq) 
errors, since the kinetic term is a dimension four operator. The size of the additional errors is estimated to be 
0{a) = 15—12% and 0{aasmq) = 22—6%. Again we expect that the dominant part of this error is eliminated by the 
tadpole improvement. 

In this work we adopt an Og-linear extrapolation for the continuum limit, because the leading order scaling violation 
is always 0((agAQCD)^, OgAQCDWg) irrespective of the choice of C. We also perform an Og-linear extrapolation to 
estimate systematic errors. In practice we use results at three finest lattice spacings i.e. f3 = 5.90-6.35 [usmq < 1) for 
the continuum extrapolation, excluding results at /3 = 5.70 (agraq > 1), which appear to have larger discretization 
errors as expected from the naive order estimate. Performing such extrapolations for all sets of A/iat = (Mpoie, Afidn) 
and ( = C^^), adopt the choice which shows the smoothest scaling behavior for the hnal value, and use others 

to estimate the systematic errors. 


IV. RESULTS 

Now we present our results of the quenched charmonium spectrum obtained with the anisotropic quark action. In 
this section, we first compare results of with Second, we determine the lattice scale, and study the effect of 
(Aflat, C) tuning. We then show the results of charmonium masses and mass splittings, and estimate their continuum 
limit. 


A. Dispersion relation and ( 

In Fig. H, we plot a typical example of the dispersion relation and the speed of light. As shown in the left figure, 
the linearity of in is satisfied well. Indeed the ‘effective speed of light’, dehned by 


Ceff(p) = 


/f;(p)2-f;(o)2 


(54) 
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has a wide plateau as shown in the right figure. Therefore we employ the linear fit in to extract the speed of light 
c from E^. This figure also illustrates that the speed of light c for agrees well with that for within errors. 
This is indeed the case for all data points as observed in Table ||. The speed of light c seems universal for all mesons 
as pointed out in Ref. Q. 

The nonperturbative value of (, is obtained by demanding that the speed of light c is equal to 1 within 1%. 
On the other hand, the tree-level tadpole-improved value, gives c deviating from 1 by 2—4% i.e. 2—4 a at most, 
which is much smaller than the size of the 0{a,aasmq) error (12—15%, 6—22%) estimated in the previous section. 
This suggests that 0{a,aasmq) errors associated with are almost eliminated by the tadpole improvement, as 
expected. 

In Fig. H, and at tuq = and TOq are plotted as a function of mo = OtmJJ. We 

find that (circles) and (squares and solid line) agree within errors at mo <0.3 but deviate from each other 
at riiQ ~ 0.5 (/3 = 5.7). The latter is one of the reasons why we exclude this point in the continuum extrapolation. 
One also notices that the slope of v approaching the value z/ = 1 in the continuum limit is steep, and in addition, the 
difference z/'^^ — for our data does not have a smooth dependence in otmj}. As discussed in Sec.^ these features 
of z/^^ bring complications in the scaling behavior of the hyperfine splitting. 


B. Lattice scale 

In this work, we determine the lattice spacing via the Sommer scale tq ]^ ], the IP — 15 meson mass splitting, and 
the 25— 15 splitting. We compare the results obtained with these different scales, in order to estimate the quenching 
errors. 


1. Scale from the Sommer scale ro 


In order to calculate the static quark potential needed for the extraction of cq, additional pure gauge simulations 
listed in Table VI are performed. Using Lag > 1.4 fm lattices, we measure the smeared Wilson loops at every 100-200 
sweeps at six values of j3 in the range /3 = 5.70—6.35. Details of the smearing method |3^j3^ are the same as those 
in Ref. |3q. We determine the potential V{f) at each (3 from a correlated fit with the ansatz 


IF(r,£) =C'(f)e“*'^(^)-‘, 


(55) 


where f = r jag and t = tjat are the spatial and temporal extent of the Wilson loop in lattice units. The fitting range 
of i is chosen by inspecting the plateau of the effective potential atVes{f,i) = log(IU(r,£)/IU(f,£ -|- 1)). A correlated 
fit to V(f) is then performed with the ansatz 


atV{f) = atVo + {atasa)f - {e/C)\ + atSV, atSV = I ( ^ - 

r \r 


where a is the string tension and [4] is the lattice Coulomb term from one-gluon exchange 

d^k cos(k • f) 


= 47r 


U (2^)3 4ELiSin"(fc*««/2). 

We extract ro/ag from the condition that |r=ro = c, i.e. 


(56) 


(57) 


ro/ag = 


c — e 

^atUga 


(58) 


with c = 1.65. The error of ro/ag is estimated by adding the jack-knife error with bin size 5 and the variation over the 
fitting range of f. Keeping to the ansatz Eq. (|^), we attempt three different fits: (i) 2-parameter fit with e = 7r/12 
and I = 0 fixed, (ii) 3-parameter fit with e = 7r/12 fixed, and (iii) 4-parameter fit. We check that ro/ag from these 
three fits agree well within errors (See Fig. 0). We adopt ro/ag from the 2-parameter fit as our final value. Results 
of ro /Og at each (3 are summarized in Table VI. 

Next we fit tq /og as a function of (3 with the ansatz proposed by Allton , 
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{aJro){P) = f{P) {1 + C2 a{Pf + C4 a{P)* )/co , 


(59) 


x{/3) = 


where (3i = 6.00 and /(/3) is the two-loop scaling function of SU(3) gauge theory, 


fiPi) ’ 


1 11 109 

m = 6/.^) ^ ^ exp(-^) , ^0 = ^ , 54 = ^ , 

and Cra(n = 2,4) parameterize deviations from the two-loop scaling. From this fit, we obtain that 

Co = 0.01230(29), C 2 = 0.163(54), C 4 = 0.053(22) 


(60) 


(61) 


with /^DF — 0.51. As shown in Fig. |T^, the fit curves reproduce the data very well. We use Eq. (0) in our later 
analysis. Finally, we obtain from the input of rp = 0.50 fm. The values of Og at each (3 are given in Table 0. 


2. Scale from charmonium mass splittings 


The quarkonium IP — 15" and 2S — 15 splittings are often used to set the scale in heavy quark simulations since 
the experimental values are well determined and they are roughly independent of quark mass for charm and bottom. 
Here we take the spin average for 15, IP and 25 masses, so that the most of the uncertainties from the spin splitting 
cancel out. The lattice spacing at mg = is given by 


«? = ^Qlat/Qe 


(Q = AM(1P - 15), AM(25 - 15)), 


(62) 


where Qiat denotes the value in the tempo ral l attice unit. We use the data of {MpoieX"^^) and check that other 
choices do not change af sizably. In Table VH we summarize the values of and aj for all Q including rp, 

and plot the /3-dependence of aj in Fig. H, We observe that < a’^° < holds for (3 = 5.70—6.35. To 

show this explicitly, on the right we also plot the ratio /dfS and / 0 ^° as a function of 0 ^°. Deviations 

at 


from unity are about —5% for -|-(10—15)% for /af° and hence -|-(10—25)% for a 


,2S-1S/„1P-1S 




our simulation points. The major source of discrepancy among the lattice spacings from different observables is the 
quenching effect. Another source is the uncertainty of input value of rp = 0.50 fm, which is only a phenomenological 
estimate. Other systematic errors are expected for for the following reasons. Our fitting for 25 masses may 

be contaminated by higher excited states. In addition, the lattice size ~ 1.6 fm may be too small to avoid finite size 
effects for 25 masses. On the other hand, the fitting for IP masses are more reliable, and we have checked that the 
finite size effects are negligible for AM(1P— 15) in preparatory simulations (see also Ref. |Q). For these reasons, we 
consider the scale to be the best choice for physical results on the spectrum. We present the results for three 

scales in the following, however, to show the dependence of the spectrum on the choice of the input for the lattice 
spacing. In order to make a comparison with the results by Klassen and Chen, who employ rp to set the scale, we 
use the results with a'y. 


C. Effect of (Miat,C) tuning 


In Fig. 12, we plot the results of spin-averaged mass splittings and spin mass splittings for each choice of (Miat, C)- 
The upper two figures show the spin-averaged splittings AM(1P— 15) and AM{2S — 15), while the lower two show 
the S-state hyperfine splitting AM(1^5i — 1 ^5o) and the P-state fine structure AM(l^Pi — l^Pp). Numerical values 
for each choice at /3 = 6.1 are given in Table Vlll , Here we set the scale with rp because it has the smallest statistical 
error. 

For all of mass splittings in Fig. the results for (Mpoie, C'^^) — (Afkin, C^^) '^ell agree with those for (Mkin, 
suggesting that the mass splittings are independent of the choice of C whenever the Mkin tuning is adopted. This 
can be understood as follows . Setting the measured kinetic mass to the experimental value Mkin = Mexp for the 
meson roughly corresponds to setting m 2 = mcharm for the quark, where the kinetic mass for the quark m 2 is given 
by Eq. ( 0 ) at the tree level. Since the spin-averaged splitting is dominated by m 2 , setting m 2 = mcharm for each C 
results in the same value for this splitting. With our choice of the spatial clover coefficient Cg = rg, ms = m ,2 also 
holds independent of C, at the tree level. Hence the ^in splitting takes approximately the same value because it is 
dominated by the magnetic mass rriB given by Eq. (|l4|). 
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As a result, we practically have only two choices for (Miat,C)i *-e- (Afpoie, and (Mpoie,C'^^) — (Mkin,C'^^) — 
(Afkin, As observed in Fig. |I^, however, the results for (MpoiejC"'"^) agree with those for the other choices at 

three finest Og, within a few a for the hyperfine splitting and lu for other mass splittings. This shows that the choice 
(Mpoie, is as acceptable as any other, with our numerical accuracy, for the lattices we adopted. Since the hyperfine 
splitting for the choice (Mnnisi C^) has a smoother lattice spacing dependence (at /3 > 5.9) and a smaller error than 
that for other choices in Fig. ^ we decide to use the data with (Mpoie, for the continuum extrapolations. The 
results for other choices are used to estimate the systematic errors. A slight bump in the lattice spacing dependence 
of the hyperfine splitting for (Mpoie, is in part ascribed to the statistical error of itself, as discussed in Sec.^. 


D. The charmonium spectrum 


The results for charmonium spectrum, obtained for (Mpoie, for the three choices of scale are plotted in Fig. 
together with the experimental values, and numerical values are listed in Tables IX-Xl. As observed in Fig. |^, the 
gross features of the mass spectrum are consistent with the experiment. For example, the splittings among the Xc 
states are resolved well and with the correct ordering (xco < Xa < Xc2)- Statistical errors for the IS, IP and 2S state 
masses are of 1 MeV, 10 MeV and 30 MeV, respectively. When we set the scale from the IP — 1A (2S — 15) splitting, 
the spin structure and the 2S — 15 (IP — l5) splittings are predictions from our simulations. 


E. S-state hyperfine splitting 


We now discuss our results for the S-state hyperfine splitting AM(l^Si — 1^5o), which is the most interesting 
quantity in this work. The hyperfine splitting (HFS), arising from the spin-spin interaction between quarks, is very 
sensitive to the choice of the clover term, as noticed from Eqs. (0) and (HI)- Since the clover term also controls the 
lattice discretization error of the fermion sector, the calculation of the HFS is a good testing ground for the lattice 
quark action. 

In Fig. 1^ we plot our results for the S-state HFS with (Mpoie, C"*"^) for each scale input by filled symbols. From the 
Og-linear continuum extrapolation using 3 points at /3 = 5.90-6.35, we obtain 


AM(13S'i - 1 ^ 50 ) = 


72.6(0.9)(-hl.2)(-3.8) MeV 
85.3(4.4)(-p5.7)(-2.5) MeV 
53.9(5.8)(-1.5)(-2.0) MeV 
117.1(1.8) MeV 


(ro_input) 

(IP — 15 input) 
(25 — l5 input) 

(experiment) 


(63) 


where the first error is the statistical error. The second error represents the ambiguity in the continuum extrapolation, 
estimated as the difference between the a^-linear and the Os-linear fits. The third error is the systematic error 
associated with the choice of (Miat, C)- We estimate it from the maximum difference at the continuum limit between 
the choice of (Mpoie, C"^^) the other three choices. Our estimate of the S-state HFS is smaller than the experimental 
value by 27 % if the IP — l5 splitting is used to set the scale. A probable source for this large deviation is quenching 
effects. 

In this figure, we also plot previous anisotropic results by Klassen (set B in Table IH) Q and Chen (set C) at 
^ = 2 and 3 with the same choice of the clover coefficients Cs,t and using rg to set the scale. The difference between 
our simulation and theirs is the choice of C and the tadpole factor for Cs,t, as noted in Table We use and the 
tadpole factor estimated from the plaquette , while they used and tadpole estimate from the mean link in the 
Landau gauge u^. As shown in this figure, our result in the continuum limit with rg input agrees with the results by 
Klassen and Chen The results with a different choice of the clover coefficients Cs,t by Klassen (set D) will 
be shown in Sec.^, where we will study the effect of Cg to the HFS. 


F. P-state fine structure 

Results for the P-state fine structure are shown in Figs. and |^. The value of the P-state fine structure in the 
continuum limit and the systematic errors are estimated in a similar manner to the case of the S-state HFS. For 
l^Pi — l^Pg splitting, we obtain 
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(64) 


AM{l^Pi - l^Po) 


68.4(5.0)(+11.8)(-3.0) MeV 
79.2(6.6)(+16.5)(-2.4) MeV 
50.5(6.2)(+7.9)(-2.2) MeV 
95.5(0.8) MeV 


(ro_ input) 

(IP — 15' input) 
(25 — 15 input), 
(experiment) 


Note that the systematic errors from the choice of the fit ansatz (second error) are rather large here, due to the 
large scaling violation seen in Fig. |^. The result with the IP — 15 input yields a 17% (2.5a) smaller value than the 
experiment. Our result with the tq input is consistent with the previous results by Klassen |19| and Chen [p4[ . 

For 1^P2 — l^Pi splitting, we obtain 


AM(1^P2 - l^Pi) = 


31.1(8.4)(+8.1)(-1.0) MeV 
35.0(9.0)(+9.6)(-0.7) MeV 
23.7(6.1)(+5.6)(-0.8) MeV 
45.7(0.2) MeV 


(ro_ input) 

(IP — 15 input) 
(25 — 15 input) 
(experiment) 


(65) 


where we use the result from the E representation operator for ^P 2 . As observed in Tables IX-XI, the mass difference 


AM(\^P 2 t — 1 ^P 2 £:) is always consistent with zero, suggesting that the rotational invariance for this quantity is 
restored well in our approach. The value of AM(1^P2 — l^Pi) is smaller than the experimental one by 23% (Icr) with 
the IP — l5 input. There is no lattice result from the anisotropic relativistic approach to be compared with. 

Next we consider the ratio of the two fine structures, AM(1^P2 — l^Pi)/AM(l^Pi — l^Po). In Fig. 0 we plot the 
lattice spacing dependence of this ratio. As shown in this figure, the scaling violation of the ratio is smaller than that 
for the individual splittings (Figs. and |^. Moreover, results are always consistent with the experimental value 
within errors. Presumably this is in part due to a cancellation of systematic errors such as the discretization effect 
and the quenching effect in the ratio. Our continuum estimate of this ratio is 


AM(13P2 - l3Pi) 
AM(l3Pi - l3Po) 


0.47(14)(+06) 

0.45(14)(+05) 

0.49(13)(+06) 

0.48(00) 


(ro_ input) 

(IP — 15 input) 
(25 — 15 input) 
(experiment) 


( 66 ) 


Our results agrees well with the experimental value. We omit the systematic error arising from the choice of (Mjat, C)j 
which is found to be much smaller than others. 

Another interesting quantity is the P-state hyperfine splitting, AM{l^Pi — l^P), where M(1^P) = [5M(l3p2) + 
3M(l^Pi) + M(l3Pg)]/9. This should be much smaller than the S-state hyperfin e sp litting because the P-state 
wavefunction vanishes at the origin. The lattice spacing dependence is shown in Fig. 18 and the continuum estimate 
is 


AM{l^Pi - l^P) 


-1.4(4.0)(-f0.6) MeV 
-1.5(4.6)(-f0.7) MeV 
-1.5(2.6)(-f0.3) MeV 
+0.9(0.3) MeV 


(ro_ input) 

(IP — 15 input) 
(25 — 15 input) 
(experiment) 


(67) 


The sign is always negative at finite Qs and in the continuum limit, but within errors the continuum value is consistent 
with the experimental value. We do not observe sizable differences between results using different scale inputs for this 
quantity. 


G. IP — IS splitting 

The mass splittings between the orbital (radial) exited state and the ground state such as the IP — 15 (25 — 15) 
splitting are dominated by the kinetic term in the non-relativistic Hamiltonian Eg. &■ Since the dependence on 
the choice of (Miat,C) is small compared to the statistical error, as shown in Eig. we ignore the systematic error 
from the choice of (Miat,C) in this and next subsections. Results of the spin-averaged and spin-dependent IP — 15 
splittings are shown in Figs. and In the continuum limit, the spin-averaged IP — 15 splitting is 

( 413(14)(—15) MeV (ro input) 

AM(1P-15) = < 351(21)(-20) MeV (25 - 15 input) . (68) 

458(01) MeV (experiment) 

The spin-dependent IP — 15 splitting deviates from the experimental value by 0-10% (1-5 (t) with the rg input and 
15-25% (3-5cr) with the 25 — 15 input, as shown in Fig. The result of the l^Pi — 15 splitting with the rg input 
agrees with the result by Chen within a few a in the continuum limit. 
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H. 2S — IS' and 2P — IP splittings 


In Figs. E and 1^, we show the results of the spin-averaged and spin-dependent 2S— IS splittings. In the continuum 
limit, these splittings deviate from the experimental values by ^20% (2.5cr) with the tq input and ~30% (4cr) with 
the IP — IS input. For the spin-averaged 2S — IS splitting, we obtain 

( 701(40)(-fl3) MeV (ro_inpu9 

AM(2S-1S) = 772(47)(-f35) MeV (IP - IS input) . (69) 

I 595(01) MeV (experiment) 


Besides quenching effects, possible sources of the deviations are finite size effects and the mixing of the 2S with higher 
excited states. Figure 23 shows the result for 2P — IP splittings. Note that there is no experimental value for this 
splitting at present. Our results of 2S — IS and 2P — IP splittings_are consistent with previous results by Chen. We 
also calculate mass splittings such as AM(2^Si — 2^So) and AM(2P — 2S), but these suffer from large statistical and 
systematic errors. We leave accurate determinations of the excited state masses for fnture studies. 


I. The charmonium spectrum in the continuum limit 


We summarize the continuum results for the charmonium spectra obtained with the data of (Mpoie, the 

Us-linear fit ansatz in Fig. where the scale is set by IP — 15 splitting. Numerical values for three scales are 
listed in Tables where the errors are only statistical. Among three different scales, results with the IP — 15 

input are the closest to the experimental value for the ground state masses. The spin splittings such as the hyperfine 
splitting AM(1^5i — 1^5o) and the fine structure AM(l^Pi — I^Pq) are always smaller than the experimental values 
irrespective of the choice of the scale input, which is considered to be quenching effects. 


V. EFFECT OF CLOVER COEFFICIENT FOR HYPERFINE SPLITTING 


We now come back to the issue of the hyperfine splitting. In Sec.tVE, we have shown that our result of the HFS 
(set A in Table ID) agrees with previous results by Klassen (set B) and Chen (set C) in the continuum limit, with the 
same choice of the clover coefficients Eqs. (^) and (^). However, as mentioned in the Introduction, when Klassen 
made a different choice of the clover coefficients (set D), he obtained apparently different values of the HFS in the 
continuum limit. This choice is given by ^ Cs = Ijv where the tilde denotes the tadpole improvement, Cg = u^gCs- 
Since v ^ 1 as asUiq —> 0, it agrees with the correct choice Cs = 1 in the limit Og ^ 0 with fixed m^, but is incorrect 
at finite Og. The quark action then generates an additional 0(asAgcDW^) error. Even with such a choice, if Ugniq is 
small enough, the result should converge to a universal value after the continuum extrapolation. However, in Refs. 
Klassen obtained HES(ag = 0, cq input) « 95 MeV with Cg = l/i^, which is much larger than the result 
= 0, To input) ~ 75 MeV with c'g = 1 both by Klassen and in the present work. 

A possible source of this discrepancy is a large mass-dependent error of 0(agAQCD • (os^g)"') (n = 1,2,- ■ •) for the 
results with Cs = l/i'. In fact, Klassen adopted rather coarse lattices with agmq « 1—2, for which such errors may 
not be negligible. Because the HES is sensitive to the spatial clover term, the choice of Cg = 1/:^ may then result in a 
non-linear Og dependence for the HFS. In the following, in order to study the effect of the choice of the spatial clover 
coefficient Cg to the HFS, we make a leading order analysis motivated by the potential model and compare it with 
numerical results, which will give us a better understanding of the above problem of the HFS. 

The potential model predicts that, at the leading order in both a and l/iriq, 


HFScont f —) • f —) l'f'(0)lco„t , (70) 

\mqj \mqj 

where ruq = uiq for the quarkonium, Sq^q are quark and anti-quark spins, and 4^(0) is the wavefunction at origin. 
HFScont is the hyperfine splitting in the continuum quenched (n/ = 0) theory, which is not necessarily equal to 
the experimental value. In non-relativistic QCD, the • Sg interaction arises from the S • B term for quark and 


^ This choice corresponds to d; = 1 in the mass form notation Eq. (^), while the correct choice 


1 corresponds to ui = u. 
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anti-quark. Giving the non-relativistic interpretation to our anisotropic lattice action, we expect that the lattice HFS 
is effectively given by 

(71) 

V ms / yruBj 

where tob is the magnetic mass Eq. (^) in the effective Hamiltonian. Therefore, in our approach, HFS is dominated 
by the magnitude of l/m^, which depends on the spatial clover coefficient Cg. The ratio, 


HFSiat |^(0)IL 

HFSeont ^ l«'(0)|2„„t 


generally deviates from 1 at finite Ug, and should approach 1 as Og ^ 0. At the leading order in a, |'I'(0)|^oat 
while oc m 2 with m 2 the kinetic mass Eq. (p^). Since m 2 does not depend on the spatial clover coefficient Cg 

at the tree level, we neglect the lattice artifact for and set |^'(0)|fat/l'f'(0)|cont = 1 in the following, which is 

sufficient for the present purpose. Now we define 


7?hfs = 



/ Of mg y 
\atmB ) 


(73) 


as a measure of lattice artifacts for the HFS, where the tilde denotes the tadpole improvement. In the continuum 
limit, i?HFS = 1- Since m^ is constant independent of Og, we identify m^ with rdi for the pole mass tuning (i.e. when 
setting the measured pole mass to the experimental value Mpoie = Afexp for the meson), and with m 2 for the kinetic 
mass tuning (Mkin = Mexp). 

At the tree level with the tadpole improvement, the pole mass mi, the kinetic mass m 2 and the magnetic mass ms 
for the quark are given by 


atrfii = 

= log(I -f mo), 


(74) 

I 

2z/2 


(75) 

atm2 

mo(2 -1- mo) 

l + rrio’ 

1 


^CsV 

(76) 

atruB 

mo(2 -1- mo) 

l + mp' 


where v = ^o/Cj Cg = itgCg, and rdo = atvigo is given by Eq. (|^). To obtain Eqs. (|7^ ) and (f^), we use the formula 
^ = ^0 = iut/us)^o- In the following we present the Cgm^ dependence of i?HFS in the case of Cg = I (set A,B,C) and 
Cg = l/v (set D), and compare them with the corresponding numerical data for the S-state HFS. For the definition 
of C (or v), there are two choices adopted so far, the tree level tadpole improved value and nonperturbative one 
At ^ rfii = W 2 for the quark, but Mpoie ^ Mkin for the measured meson. On the other hand, at 

mi ^ m 2 though Mpoie = Afkin- Thus in the case of ^ ^poie = Mkin tuning, the identification of m^ 

(=mi or m 2 ) in i?HFS Eq- (^H) mentioned above is ambiguous. Although such an ambiguity should vanish in the 
continuum limit, we present i?HFS with both m^ = mi and m^ = m 2 to check consistency. For actual numerical data 
of the HFS, we focus on the results with the ro input because Klassen has adopted rg for the scale setting. 


A. The case oi Cs = l/v 


First we consider the case of Cs = Ijv (set D), which is correct only for Ogm^ = 0 at the tree level. In Fig. we 
plot the (osmg)^ dependence of i?HFS at ^ = 3 and 2 for c'g = Ijv with v = Numerical values of 

were taken from Ref. |I9|. Because of the ambiguity for mq mentioned above, we show the results with m^ = rhi and 
mg = m 2 ; the difference between them decreases as 0, as expected. We have checked that plotting Rhfs as a 
function of Og, instead of (ogmg)^, does not change the figure qualitatively. We also plot the results with Cg = Ijv but 
v = = Co/C"'"^ where mi = rn 2 holds, as a dotted line (^ = 3) and a dashed line (^ = 2) for a guide to the eye. As 

shown in this figure, Rhfs has a non-linear aj. dependence toward the continuum limit (=1), indicating that the mass 
dependent error is large for the region Ogmg = 1—2. Rhfs is larger than I even at (ogmg)^ ~ I, which suggests that 
the actual HFS should rapidly decrease toward Og ^ 0, and data at (ogmg)^ < I are needed for a reliable continuum 
extrapolation for the HFS. 

Now let us compare Rhfs with numerical results of HFS. In Fig. we plot corresponding results of HFS by 
Klassen for Cg = I/z^ 0- The results at ^ = 3 for Cg = l/v are clearly larger than the results for c'g = I (see filled 
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circles in Fig. 
.,2 


and the results at ^ = 3 and 2 


appear to converge to « 95 MeV in the continuum limit with 
we find that the lattice spacing dependence of the 


an a“-linear scaling. However, comparing Fig. ^ and Fig. 
numerical data of HFS qualitatively agrees with that of i?HFs: for both HFS and i?HFS; data at ^ = 3 are larger than 
data at ^ = 2, and the difference between ^ = 3 and 2 decreases as Os —> 0. From an a^-linear extrapolation of i?HFS 
using the finest three data points, we obtain i?HFS ~ 1.2—1.3 at = 0. Because the correct continuum limit of i?HFS 
is 1, this suggests a 20—30% overestimate from the neglect of non-linear dependence of i?HFS on a^. Hence the result 
with c's = l/j 2 , HFS(as = 0) ~ 95 MeV, reported in Refs. |l^J^ is likely an overestimate by 20—30%. 

These analyses indicate that the origins of this overestimate are, first, the choice for the spatial clover coefficient 
c's = l/u (= l/i/^^), and second, the use of coarse lattices with agniq > 1. As shown in Fig. ||, i/ (= 1/cs in this case) 
should eventually start to move up to 1 linearly around atmj}<0.3, which corresponds to (asmg)^<0.6 in Fig. but 
Klassen’s data of (open diamonds) do not reach such a region. We conclude that the continuum extrapolation 
for the HFS should not be performed using the data on such coarse lattices, and results at finer lattice spacing are 
required. 


B. The case of ct, = 1 


Next we consider the case of Cg = 1 (set A, B and C), whic h is c orrect for any asTUq at the tree level. In this case, 
there are two choices for C, and As mentioned in Sec. IV C, tub = 'hi 2 holds for both choices of with c's = 1. 

In the case of C = C'^^ which is adopted only in our work (set A) so far, Rrfs = I is always satisfied, since 
nil = = rn-B by definition. This suggests that the scaling violation of HFS for c's = I should be much smaller than 

that for c's = Xjv. The numerical result for the HFS with the pole mass tuning has already been shown in Fig. [l^ 
and re-plotted in Fig. by filled circles, which gives our best estimate, HFS(as = 0) = 73 MeV. 

We next consider the case of C = where Mpoie = Afkin for the measured meson. When we identify rUg = TO 2 , 
Rhfs = 1 is always satisfied again because m 2 = tob even at ^ When we identify niq = toi, Rhfs ^ 1 in 

general, due to the deviation of from The results of Rhfs with mq = mi at C = are shown in Fig. |2^, 
and corresponding numerical results for the HFS are shown in Fig. Comparing Fig. with Fig. ^we again note 
that the lattice spacing dependence of the HFS qualitatively agrees with that of RhfSi *-e-) for both HFS and Rhfs, 
data at f = 3 by Klassen (open diamonds, set B) and those at ^ = 2 by Chen (open triangles, set C) are close to each 
other and larger than our data at C = ag-linear extrapolation using the finest three data points gives HFS 

~ 70—75 MeV and Rhfs ~ 0.9—1.0 at Ug = 0. The latter confirms that a continuum estimate of HFS with Cg = I is 
more reliable than that with Cg = l/v. 

Concerning our results at ^ = 3, as shown in Fig. Rhfs for C = (stars) does not scale smoothly around 
(agmg)^<l, while that for C = (filled circles) is always unity. This behavior is caused by the fact that the difference, 
^NP _ monotonic in Ogm, (see Fig. |^). Correspondingly the numerical value of the HFS, displayed in Fig. 

also shows a slightly non-smooth lattice spacing dependence near ~ 0, which qualitatively agrees with the (ogm^R 
dependence of Rhfs in this region. A possible source of this behavior is the statistical error of itself, because 
HFS (i?HFs) is also sensitive to the value of C as well as Cg. Due to this reason, we have not used the results with 
^ main analysis in Sec. IV. 


VI. CONCLUSION 

In this article, we have investigated the properties of anisotropic lattice QCD for heavy quarks by studying the 
charmonium spectrum in detail. We performed simulations adopting lattices finer than those in the previous studies 
by Klassen and Chen, and made a more careful analysis for 0((agmg)") errors. In addition, using derivative operators, 
we obtained the complete P-state fine structure, which has not been addressed in the previous studies. 

From the tree-level analysis for the effective Hamiltonian, we found that the mass dependent tuning of parameters 
is essentially important. In particular, with the choice of rg = 1 for the spatial Wilson coefficient, an explicit Ugm^o 
dependence remains for the parameters C, and ct even at the tree level. Moreover we have shown in the leading order 
analysis that, unless the spatial clover coefficient Cg is correctly tuned, the hyperfine splitting has a large 0((agmg)"') 
errors, which can explain a large value of the hyperfine splitting in the continuum limit from rather coarse lattices in 
the previous calculation by Klassen. On the other hand, if c~g is mass-dependently tuned, the continuum extrapolation 
is expected to be smooth for the hyperfine splitting. 

Based on these observations, we employed the anisotropic clover action with rg = 1 and tuned the parameters 
mass-dependently at the tree level combined with the tadpole improvement. We then computed the charmonium 
spectrum in the quenched approximation on ^ = 3 lattices with spatial lattice spacings of OgTOg < I. A fine resolution 
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in the temporal direction enabled a precise determination of the masses of S- and P- states which is accurate enough 
to be compared with the experimental values. Our results are consistent with previous results at ^ = 2 obtained by 
Chen 1^^, and the scaling behavior of the hyperfine splitting is well explained by the theoretical analysis. We then 
conclude that the anisotropic clover action with the mass-dependent parameters at the tadpole-improved tree level is 
sufficiently accurate for the charm quark to avoid large discretization errors due to heavy quark. We note, however, 
that agtriq < 1 is still necessary for a reliable continuum extrapolation. 

We found in our results that the gross features of the spectrum are consistent with the experiment. Quantitatively, 
however, the S-state hyperfine splitting deviates from the experimental value by about 30% (7cr), and the P-state fine 
structure differs by about 20% (2.5cr), if the scale is set from the IP — 15 splitting. We consider that a major source 
for these deviations is the quenched approximation. 

Certainly further investigations are necessary to conclude that the anisotropic QCD can be used for quarks heavier 
than the charm. In particular it is important to determine the clover coefficients as well as other parameters non- 
perturbatively, since the spin splittings are very sensitive to the clover coefficients. It is also interesting to calculate 
the spectrum with Vg = 1/^ and compare the result with the current one in this paper, since the notorious agmqo 
dependence vanishes from the parameters with this choice at the tree level. Finally full QCD calculations includ¬ 
ing dynamical quarks are needed to establish the theoretical prediction without systematic errors for an ultimate 
comparison with the experimental spectrum. 
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APPENDIX A: DERIVATION OF HAMILTONIAN ON THE ANISOTROPIC LATTICE 


The lattice Hamiltonian H is identified with the logarithm of the transfer matrix T: 

H = -logf. 


(Al) 


T and H for the asymmetric clover quark action on the isotropic lattice have been derived in Ref. [pd] . An extension to 

the anisotropic lattice is straightforward. Using the fields 'k and ^ = 'k^7o which satisfy canonical anti-commutation 
relations, the Hamiltonian in temporal lattice units H for the anisotropic quark action is given by 


= 4' 


0(7711 - 


Cpa^ 


-(r^D^ + iCs'S ■ B) - <F/i(TOo)as0 - CF/2(wo)a^0^ 


2(1 -I- TTlo) 

where {T,i,ai) = {-^e^jkO-jk,-icroi), {Bi,Ei) = {\eijkFjk,Foi) and 

0(7771 = log(l -I- Too), 


4' + 0(p"o^), 


and 


0 = 7(7 • D -k -(1 - C()o(a • E), 


4 ^ 2(l-kx)log(l-kx) ^ fl{x) 

mx) = - 77 ——^-, f 2 (x) = 


x(2 + x) ’ 2 log(l + x) x(2 -k x) 

Therefore the lattice Hamiltonian in physical units is given by 


1 

at 


-H = ^ 


= ^ 


mi — 


mi — 


2(1 -k ?77o) 


(r^D^ -k iCsT, ■ B) - 7Cf/i(wo)Co0 - Cy2{mo)CoatO" 


4- + 0(p"af) 


-(r'D2 + 7c;S.B)-7C^/i(mo)0-CF"/2(mo)o(02 ^-kO(p3o2), 


2(1 -k mo) 


(A2) 

(A3) 

(A4) 

(A5) 

(A6) 

(A7) 
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where 


Cf = ■CoCf, r' = ^oTs, c'^ = ^qCs. 


(A 8 ) 


Note that Eq. (A7) for the anisotropic lattice is the same as that for the isotropic lattice except for use of {ot, r', Cg} 
instead of {a, Cfj ^s, Cs}. Thus one can repeat the derivation of the tree level value of bare parameters (Cf and Cs,t) in 
Ref. even for the anisotropic lattice, after replacing {a, Cfj Cg} by {ot, Cp, c(,}. 

When the lattice Hamiltonian is expressed in more continuum-like form 


1 


at 


H = \It[ 6 oTOq -|- 6 i7 ■ D -|- 04620 ^ -f iatbs^ ■ B -|- atbECt ■ E -b atbsojoll • D ,7 • E] -I- 


(A9) 


the coefficients b are given by 


bo = milmq, (AlO) 

= CF/i("io), (All) 

‘^ = -2n^ + AA(mo), (A12) 

bs = ^(1 - Ct)CF/i(mo), (A14) 

&SO =-^(1 - Ct)CF^/ 2 (nio)- (A15) 


In order to determine tree level parameters, the lattice Hamiltonian should be matched to the continuum one to the 
desired order in Ug. The continuum Hamiltonian to which the lattice one is matched is either the Dirac Hamiltonian 
l?Dirac = at^{niq -b 7 • D)'k, Or the non-relativistic Hamiltonian Hnr = at^{mq + 70 Aq — -b • • •)'!'• Both choices 
give the same tree level parameters. 

In the Hamiltonian formalism, the unitary transformation U is possible because the eigenvalues of H are invariant 
under it. For example, consider a unitary transformation, 

# ^ 1/4', ^ (A16) 


with 


U = exp(—at 0 i 7 • D — a^OEOt ■ E), 


(A17) 


where 0i and Oe are parameters. This is called the Foldy-Wouthuysen-Tani (FWT) transformation, whose element is 
a spin off-diagonal matrix. After this transformation the coefficients b become 


b^ = bo, (A18) 

b^ = bi- 2mqatho9i, (A19) 

6 ^ = 62 - 2bi6i -b 2mqatboel, (A20) 

b% = bB - 2bi9i + 2mqatbo9l, (A21) 

b% = bE - 9i - 2mqatbo9E, (A22) 

b^o = bso - ^9^ + bE9i + bi9E - 2mqatbo9i9E. (A23) 


The transformed Hamiltonian with b^ is matched to either i7Dirac or TInr so as to obtain tree level parameters. 
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P 

C 

^0 

Cs 

Ct 

a? [fm] 


X 

T 

Las [fm] 

5.70 

3 

2.346 

1.966 

2.505 

0.204 

8'^ 

X 

48 

1.63 

5.90 

3 

2.411 

1.840 

2.451 

0.137 

12® 

X 

72 

1.65 

6.10 

3 

2.461 

1.762 

2.416 

0.099 

16® 

X 

96 

1.59 

6.35 

3 

2.510 

1.690 

2.382 

0.070 

24® 

X 

144 

1.67 


TABLE 1. Simulation parameters. Las is calculated using a^”, the lattice spacing determined from ro. 
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p 


X T 

atrriqo 

c 

sweep / conf 

#conf 

cps 

Cv 

5.70 

S'" 

X 

48 

0.320 

2.88 (NP) 

100 

1000 

1.005(10) 

1.008(11) 

5.70 

8" 

X 

48 

0.253 

2.85 (NP) 

100 

1000 

1.005(10) 

1.008(11) 

5.70 

8" 

X 

48 

0.320 

3.08 (TI) 

100 

1000 

0.962(9) 

0.965(10) 

5.70 

8" 

X 

48 

0.253 

3.03 (TI) 

100 

1000 

0.966(9) 

0.969(10) 

5.90 

12® 

X 

72 

0.144 

2.99 (NP/TI) 

100 

1000 

0.991(8) 

0.993(9) 

5.90 

12® 

X 

72 

0.090 

2.93 (NP/TI) 

100 

1000 

0.991(8) 

0.994(9) 

6.10 

16® 

X 

96 

0.056 

3.01 (NP) 

200 

600 

0.997(9) 

0.997(9) 

6.10 

16® 

X 

96 

0.024 

2.96 (NP) 

200 

600 

0.997(9) 

0.997(9) 

6.10 

16® 

X 

96 

0.056 

2.92 (TI) 

200 

600 

1.017(9) 

1.018(9) 

6.10 

16® 

X 

96 

0.024 

2.88 (TI) 

200 

600 

1.017(9) 

1.016(10) 

6.35 

24® 

X 

144 

-0.005 

2.87 (NP/TI) 

400 

400 

1.006(11) 

1.011(11) 

6.35 

24® 

X 

144 

-0.035 

2.81 (NP/TI) 

400 

400 

1.007(12) 

1.009(11) 


TABLE II. Simulation parameters continued. In fourth column, ‘NP’ and ‘TP denote the nonperturbative and tree level 
tadpole improved values for ( respectively. cps,v are the speed of light obtained from the fit for the pseudoscalar (^S'o) and 
vector (^Si) mesons. 


set 



c 

Cs 

Ct 


'^s,t 

Mi,t(15) 

scale input 

HPS(a, = 0,ro) 

A. this work 

3 

TI(m > 0), NP 

TI(m > 0) 

TI(m = 

0) 


Mpole, Mkin 

ro, IP - 15, 25 - 15 

« 75 MeV 

B. Klassen 


2,3 

NP 

TI(m > 0) 

TI(m = 

0) 


Mpolei— ALkin) 

ro 

« 75 MeV 

C. Chen |4 


2 

NP 

TI(m > 0) 

TI(m = 

0) 


Afpole(— ATkin) 

ro 

Ri 75 MeV 

D. Klassen 


2,3 

NP 

TI(m = 0) 

TI(m = 

0) 


ALpole(— ALkin) 

ro 

« 95 MeV 


TABLE III. Comparison of simulation parameters in various anisotropic lattice studies of the cc spectrum. In the third to 
fifth columns, TI(m > 0), TI(m = 0) and NP respectively denote the tree level tadpole improved value for massive quarks, that 
are correct only in the massless limit and the non-perturbative value. The sixth column shows which method is used for the 
estimation of the tadpole factors Us,t (the plaquette prescription or the Landau mean link prescription u^). The seventh 
column shows which 1,5 mass is tuned to the experimental value. The eighth column denotes quantities used for the scale 
setting. The final column is the continuum estimate of the hyperfine splitting from the a^-linear fit with the scale set by ro. 


2S+l^^ 


name 

r operator 

PA operator 


0-+ 

Vc 

^/)757/) 


®5i 

1~~ 

J/V> 

Ipji'tp 


^Pi 

1+- 

he 

ipaiji; 


"Po 

0++ 

XcO 

iptp 


®Pi 

1++ 

Xcl 

ipXi'ysi’ 


®P2 

2++ 

Xc2 


ipi-jiAi - 7jAj}V> (E rep) 
-ipi-ytAj + 'yjAi}ip (T rep) 


TABLE IV. S- and P-state operators. In the first and second columns, the state is labeled by and respectively. 

The third column shows the particle name for the charmonium family. In the fourth and fifth columns, we give the corresponding 
P operator and PA operator. 
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state 

fit form 

source 


fit range (tmis 

i/^max) 





P = 5.70 

P = 5.90 

P = 6.10 

P = 6.35 

IS,2S 

2-cosh 

00+01+11 

11/24 

17/36 

22/48 

32/72 

1P,2P 

2-cosh 

00+11+02+12 

7/18 

11/25 

15/35 

21/50 

IS, AS 

1-cosh 

01 

13/24 

19/36 

26/48 

38/72 

lS(p / 0) 

1-cosh 

01 

13/22 

20/32 

26/45 

40/66 

AP 

1-cosh 

12 

11/18 

17/25 

23/35 

33/50 


TABLE V. Fit ranges we adopted. In the first column, AS and AP denote the S- and P-state spin mass splitting respectively. 


P 

ro/tts 

LP xT 

Las [fm] 

smear# 

conf 

sweep/conf 

5.70 

2.449(35) 

12'' X 72 

2.45 

4 

150 

100 

5.90 

3.644(36) 

12® X 36 

1.65 

5 

220 


6.00 

4.359(51) 

12® X 48 

1.38 

6 

150 

100 

6.10 

5.028(35) 

16® X 48 

1.59 

6 

150 

100 

6.20 

5.822(33) 

16® X 64 

1.37 

10 

220 

100 

6.35 

7.198(52) 

24® X 72 

1.67 

12 

150 

200 


TABLE VI. Simulation parameters and results for the Sommer scale ro. The hfth column shows the number of smearing 
steps we adopted. 


P 

fo 


IP- 

IS 

2S- 

IS 



[fm] 

mg*'”"' 


mg*'”'" 

af-i^[fm] 

5.70 

0.2843(3) 

0.2037(0) 

0.2994(115) 

0.2077(30) 

0.3782(190) 

0.2272(45) 

5.90 

0.1106(2) 

0.1374(0) 

0.0972(58) 

0.1333(18) 

0.1664(150) 

0.1544(44) 

6.10 

0.0319(1) 

0.0991(0) 

0.0155(60) 

0.0934(21) 

0.0632(110) 

0.1099(37) 

6.35 

-0.0179(1) 

0.0697(0) 

-0.0301(43) 

0.0650(18) 

0.0115(84) 

0.0808(30) 


TABLE VII. Bare charm quark mass and lattice spacing a? for Q = ro, IP — IS and 2S — IS. 


(Mlat,C) 

AM(1P - IS) 

AM(2S - IS) 

AM(l®Si - CSo) 

AM(l®Pi - l®Po) 

(Mpole,^'*) 

426.7(104) 

676(30) 

71.6(07) 

57.3(37) 

(Mpo,e,C^'’) 

423.1(096) 

671(29) 

68.8(06) 

55.3(34) 

(Mkin,C'^*) 

424.1(097) 

671(31) 

69.2(14) 

55.2(38) 

(Mkin,C'"'’) 

423.6(097) 

672(30) 

69.2(13) 

55.7(37) 


TABLE VIII. Comparison of mass splittings for different choices of (Miat,C) at /I = 6.10. The results are presented in units 
of MeV, and the scale is set by ro. 
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state 

/3 = 5.70 

13 = 5.90 

f3 = 6.10 

/3 = 6.35 

as 

Exp. 


3020.9(7) 

3013.8(8) 

3014.0(10) 

3012.7(9) 

3012.7(11) 

2979.8 


3082.0(7) 

3083.1(8) 

3085.1(8) 

3083.7(8) 

3084.6(10) 

3096.9 

l^Pi 

3526.6(79) 

3506.7(57) 

3489.7(66) 

3483.8(83) 

3474.2(94) 

3526.1 

l^Po 

3496.0(94) 

3462.4(65) 

3438.7(58) 

3420.2(86) 

3408.5(95) 

3415.0 

l®Pi 

3526.7(84) 

3506.6(61) 

3490.5(62) 

3480.8(80) 

3472.3(91) 

3510.5 

1^P2E 

3555.2(106) 

3515.6(116) 

3509.8(199) 

3506.7(219) 

3503.6(250) 

3556.2 

1^P2T 

3555.0(100) 

3512.4(115) 

3508.9(179) 

3502.5(213) 

3501.2(238) 

3556.2 

IS 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6 

IP 

3536.0(85) 

3506.7(73) 

3494.0(104) 

3487.3(120) 

3480.4(137) 

3525.5 

l^Pi - IS 

459.9(79) 

440.9(59) 

422.4(67) 

417.8(84) 

407.2(95) 

458.5 

l®Po - IS 

429.2(93) 

396.7(66) 

371.3(61) 

354.2(87) 

341.2(97) 

347.4 

l®Pi - IS 

459.9(84) 

440.9(62) 

423.2(64) 

414.9(81) 

405.2(93) 

442.9 

1®P2 - IS 

488.5(106) 

449.9(117) 

442.5(198) 

440.7(218) 

436.6(249) 

488.6 

IP - IS 

469.3(85) 

441.0(74) 

426.7(104) 

421.3(121) 

413.4(138) 

457.9 

l®Si - l^So 

61.9(4) 

70.4(6) 

71.6(7) 

72.0(8) 

72.6(9) 

117.1 

l^Pl - l^Po 

32.3(34) 

46.7(34) 

57.3(37) 

62.7(42) 

68.4(50) 

95.5 

1®P2 - l®Pl 

18.1(43) 

18.2(41) 

20.4(68) 

30.4(72) 

31.1(84) 

45.7 

1®P2T — 1^P2£ 

-0.8(23) 

-2.3(28) 

-2.6(33) 

-2.0(41) 

-2.2(47) 

0.0 

l^Pl - l^P 

-6.0(18) 

-3.5(21) 

-0.7(29) 

-3.5(36) 

-1.4(40) 

0.9 

l3p, _l3p„ 

0.56(13) 

0.39(9) 

0.36(12) 

0.49(11) 

0.47(14) 

0.48 

2^% 

3719(22) 

3700(28) 

3699(32) 

3746(40) 

3739(46) 

3594 

23Si 

3767(20) 

3773(27) 

3758(31) 

3786(34) 

3777(40) 

3686 

2^Pi 

4248(68) 

4411(70) 

4214(70) 

4161(79) 

4053(95) 

- 

2^Po 

4175(93) 

4226(89) 

4148(94) 

4049(100) 

4008(122) 

- 

2^Pi 

4228(75) 

4388(77) 

4256(90) 

4140(84) 

4067(105) 

- 

2^P2e 

4238(109) 

4254(99) 

4190(144) 

4023(148) 

3992(175) 

- 

2®P2t 

4230(111) 

4281(100) 

4223(157) 

4082(146) 

4047(177) 

- 

2S 

3755(20) 

3755(27) 

3744(30) 

3776(34) 

3768(40) 

3663 

2P 

4233(74) 

4324(68) 

4209(86) 

4089(86) 

4027(105) 

- 

2P-2S 

478(73) 

569(70) 

466(90) 

313(88) 

256(107) 

- 

23Si -2iSo 

48(9) 

74(16) 

60(17) 

40(22) 

34(25) 

92 

2iSo - l^So 

698(22) 

686(28) 

685(32) 

733(40) 

726(46) 

614 

2® Si - l®Si 

685(20) 

690(27) 

673(31) 

702(34) 

692(40) 

589 

2iPi - l^Pi 

721(68) 

904(69) 

724(69) 

678(79) 

579(94) 

- 

2^Po - l®Po 

679(95) 

763(90) 

709(95) 

629(103) 

601(124) 

- 

2®Pi - l®Pi 

701(76) 

881(77) 

766(90) 

659(84) 

595(105) 

- 

2®P2 - 1®P2 

683(109) 

738(93) 

681(129) 

516(136) 

490(160) 

- 

2S- IS 

688(20) 

689(27) 

676(30) 

710(34) 

701(40) 

595 

2P- IP 

697(75) 

817(66) 

715(81) 

602(83) 

547(100) 

- 

TABLE IX. 

Results of charmonium 

masses M and 

mass splittings AM 

in units of MeV at C 

■ = using the pole mass 


tuning. The scale is set by tq. 
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state 

(3 = 5.70 

P = 5.90 

P = 6.10 

P = 6.35 

Us ^ 0 

Exp. 


3023.0(16) 

3010.3(16) 

3007.1(27) 

3004.3(33) 

3003.0(35) 

2979.8 

l^Si 

3081.4(8) 

3084.0(10) 

3087.1(12) 

3086.0(12) 

3087.5(14) 

3096.9 

l^Pi 

3515.6(29) 

3523.3(46) 

3520.7(88) 

3519.9(98) 

3518.6(106) 

3526.1 

l^Po 

3486.6(49) 

3476.2(51) 

3464.0(91) 

3446.4(92) 

3441.6(104) 

3415.0 

l^Pi 

3515.8(35) 

3523.5(44) 

3522.3(96) 

3516.8(102) 

3516.8(112) 

3510.5 

l^P2E 

3543.2(40) 

3532.9(60) 

3541.3(128) 

3544.9(139) 

3548.9(151) 

3556.2 

1^P2T 

3543.0(38) 

3529.3(69) 

3539.8(122) 

3540.0(155) 

3546.0(160) 

3556.2 

IS 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6 

IP 

3524.7(7) 

3523.4(7) 

3525.0(9) 

3523.4(8) 

3524.1(9) 

3525.5 

l^Pi - 15 

448.8(29) 

457.8(46) 

453.6(89) 

454.3(100) 

452.0(108) 

458.5 

l^Po - 15 

419.8(47) 

410.6(51) 

396.9(93) 

380.9(95) 

375.2(106) 

347.4 

l^Pi - 15 

448.9(34) 

457.9(44) 

455.3(98) 

451.3(104) 

450.3(114) 

442.9 

1®P2 - 15 

476.4(40) 

467.4(58) 

474.2(126) 

479.4(136) 

482.4(148) 

488.6 

IP - 15 

457.9(0) 

457.9(0) 

457.9(0) 

457.9(0) 

457.9(0) 

457.9 

1^51 - 1^50 

59.2(18) 

74.9(21) 

80.4(34) 

82.7(42) 

85.3(44) 

117.1 

l3pi _ l3p^ 

30.6(37) 

49.9(39) 

64.6(45) 

72.6(65) 

79.2(66) 

95.5 

l^Pz - l®Pl 

17.4(41) 

19.2(43) 

22.3(75) 

34.7(81) 

35.0(90) 

45.7 

1^P2T — 1^P2B 

-0.8(22) 

-2.5(30) 

-3.2(39) 

-2.1(51) 

-2.7(53) 

0.0 

l^Pl - 1®P 

-5.9(17) 

-3.7(22) 

-0.8(35) 

-3.7(44) 

-1.5(46) 

0.9 

l3Pi-l3Pn 

0.57(12) 

0.39(9) 

0.35(13) 

0.48(12) 

0.45(14) 

0.48 

2^% 

3704(22) 

3722(30) 

3746(39) 

3801(45) 

3806(50) 

3594 

2^51 

3749(21) 

3800(29) 

3811(41) 

3847(43) 

3849(49) 

3686 

2^Pi 

4217(70) 

4458(75) 

4294(79) 

4238(87) 

4159(100) 

- 

2^Po 

4146(95) 

4260(95) 

4222(105) 

4121(124) 

4114(138) 

- 

2^Pi 

4196(78) 

4434(83) 

4339(100) 

4222(96) 

4179(114) 

- 

2^P2e 

4203(107) 

4303(96) 

4263(145) 

4096(155) 

4091(173) 

- 

2^ P 2 T 

4194(111) 

4329(98) 

4287(163) 

4147(153) 

4131(177) 

- 

25 

3738(21) 

3781(29) 

3794(39) 

3836(42) 

3839(47) 

3663 

2P 

4200(76) 

4371(68) 

4286(81) 

4165(88) 

4132(100) 

- 

2P-2S 

462(72) 

590(72) 

492(95) 

329(97) 

290(112) 

- 

2^51 - 2^50 

45(9) 

78(18) 

65(20) 

47(27) 

43(29) 

92 

2^50 - 1^5o 

681(23) 

712(30) 

738(40) 

797(46) 

803(51) 

614 

2^51 - 1^51 

668(21) 

716(29) 

723(40) 

762(43) 

762(48) 

589 

2iPi - l^Pi 

701(69) 

935(73) 

773(76) 

718(84) 

641(97) 

- 

2^Po - l^Po 

659(96) 

783(96) 

758(106) 

674(122) 

671(137) 

- 

2^Pi - l^Pi 

681(77) 

910(82) 

817(99) 

705(94) 

662(111) 

- 

23 P 2 _ i3p2 

660(107) 

770(93) 

722(135) 

551(147) 

543(164) 

- 

25-15 

671(21) 

715(28) 

727(39) 

770(42) 

772(47) 

595 

2P - IP 

675(76) 

847(68) 

761(81) 

641(87) 

608(100) 

- 


TABLE X. 


The same as Table but the scale is set by IP 


IS splitting. 
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state 

P = 5.70 

P = 5.90 

P = 6.10 

P = 6.35 

tts ^ 0 

Exp. 


3032.3(21) 

3026.4(30) 

3024.9(33) 

3028.6(38) 

3027.4(45) 

2979.8 

l®5i 

3079.1(8) 

3079.8(10) 

3082.0(13) 

3079.5(12) 

3080.5(15) 

3096.9 

l^Pi 

3467.1(113) 

3446.7(139) 

3440.5(158) 

3415.3(170) 

3412.6(208) 

3526.1 

l^Po 

3445.3(112) 

3412.8(124) 

3398.6(130) 

3370.2(128) 

3361.5(165) 

3415.0 

l^Pi 

3467.8(117) 

3446.1(142) 

3440.1(158) 

3412.4(168) 

3409.7(207) 

3510.5 

1^P2E 

3490.4(124) 

3453.4(153) 

3460.0(198) 

3433.8(200) 

3437.7(244) 

3556.2 

1^P2T 

3490.1(120) 

3451.6(155) 

3460.0(185) 

3431.2(180) 

3435.3(226) 

3556.2 

15 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6(0) 

3067.6 

IP 

3475.2(114) 

3446.5(140) 

3445.0(164) 

3418.5(170) 

3418.2(209) 

3525.5 

l^Pi - 15 

399.7(114) 

380.2(141) 

372.8(159) 

348.5(172) 

345.1(210) 

458.5 

l^Po - 15 

377.9(113) 

346.4(126) 

330.8(131) 

303.4(131) 

294.2(168) 

347.4 

l^Pi - 15 

400.4(118) 

379.7(144) 

372.3(159) 

345.6(171) 

342.2(210) 

442.9 

1®P2 - 15 

423.0(126) 

386.9(155) 

392.2(199) 

367.0(202) 

370.4(246) 

488.6 

IP - 15 

407.8(116) 

380.1(142) 

377.3(164) 

351.7(173) 

350.8(212) 

457.9 

l®5i - 1^50 

47.4(25) 

54.4(38) 

57.7(43) 

51.5(48) 

53.9(58) 

117.1 

l3pi _ l3p^ 

23.2(29) 

35.2(35) 

45.8(46) 

43.9(54) 

50.5(62) 

95.5 

l^Pa - l^Pi 

14.1(32) 

14.4(30) 

17.3(51) 

22.2(52) 

23.7(61) 

45.7 

1^P2T — i^P2E 

-1.0(15) 

-1.7(17) 

-1.6(23) 

-1.9(24) 

-1.8(29) 

0.0 

l^Pl - l^P 

-5.4(12) 

-2.7(14) 

-0.6(21) 

-3.0(23) 

-1.5(26) 

0.9 

l3Pi-l3Pn 

0.60(12) 

0.41(8) 

0.38(11) 

0.50(10) 

0.49(13) 

0.48 

2^% 

3637(6) 

3618(8) 

3624(10) 

3641(11) 

3644(13) 

3594 

2^51 

3671(2) 

3676(3) 

3676(3) 

3669(4) 

3669(4) 

3686 

2^Pi 

4078(59) 

4241(69) 

4087(70) 

4015(76) 

3930(95) 

- 

2^Po 

4020(77) 

4103(76) 

4031(80) 

3914(88) 

3877(108) 

- 

2^Pi 

4057(66) 

4222(73) 

4125(82) 

3985(81) 

3929(103) 

- 

2^P2e 

4049(85) 

4078(85) 

4076(120) 

3884(106) 

3872(134) 

- 

2^ P 2 T 

4037(87) 

4109(84) 

4120(128) 

3958(104) 

3948(133) 

- 

25 

3663(1) 

3662(1) 

3663(1) 

3662(1) 

3663(1) 

3663 

2P 

4056(61) 

4157(65) 

4087(79) 

3945(73) 

3900(93) 

- 

2P-25 

393(61) 

495(65) 

424(79) 

283(73) 

237(93) 

- 

2^51-2150 

34(7) 

59(11) 

52(13) 

29(14) 

26(17) 

92 

2^50 - 1^5o 

605(5) 

592(8) 

600(10) 

612(10) 

616(13) 

614 

2=^5i - 1^51 

592(2) 

597(3) 

594(3) 

590(3) 

588(4) 

589 

2iPi - l^Pi 

611(57) 

794(63) 

647(64) 

600(73) 

517(88) 

- 

2^Po - l^Po 

575(77) 

690(74) 

633(79) 

543(86) 

514(105) 

- 

2^Pi - l^Pi 

589(64) 

776(67) 

685(78) 

573(76) 

520(96) 

- 

2^P2 - 1®P2 

559(85) 

624(77) 

616(109) 

450(104) 

443(128) 

- 

25-15 

595(0) 

595(0) 

595(0) 

595(0) 

595(0) 

595 

2P- IP 

581(60) 

710(58) 

642(72) 

526(70) 

487(87) 

- 


TABLE XL 


The same as Table IX 


but the scale is set by 2S 


15 splitting. 
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v(mo,^) g v(mo,^) 



FIG. 1. Tree level full mass dependences of v and ct for Va 
lits miOs = ^log(l + ?no). Vertical axis is normalized to be 



1/^ = l/^Q. Horizontal axis is the pole mass in spatial lattice 
in the massless limit. 




FIG. 2. The same as Fig. but for Va = 1. 
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FIG. 3. S-state effective masses at = 5.90, atnigo = 0.144 and ^ = 2.99. The left figure shows the l^S'o masses at p 7^ 0 
while the right shows the l^So and masses for the source ss' — 00, 01 and 11. 


1 rope. 


1 TAope. 


FIG. 4. P-state (l^A) effective masses at = 5.90, atrugo = 0.144 and = 2.99. The left figure shows the masses from the 
F operator, while the right shows those from the FA operator. 
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• l'So (1-cosh fit, 01) 

■ I^S, (1-cosh fit, 01) 

Ol'S„ (2-cosh fit, 00+01+11) 
□ I'^S, (2-cosh fit, 00+01+11) 
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FIG. 5. Fit range (tmin) dependence of masses at /3 = 5.90, atniqo = 0.144 and C, = 2.99. The legend denotes the state (fit 
ansatz, quark source). 
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FIG. 7. Bin size dependence of jack-knife error of a 
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FIG. 8. Dispersion relation (left) and speed of light (right) of S-state at /3 = 5.90, atmqo = 0.144 and C, = 2.99. On the right, 
we show the effective speed of light Ceff(p) and the speed of light c from the ht. 
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FIG. 9. The tadpole improved bare mass mo = atm^o versus = ^o/C at ^ = 3. ‘TF and ‘NP’ denote the tree level tadpole 
improved value and nonperturbative value respectively. Gircles and squares are our data at mo = mj,m§ (ss mj)*'”"') for 
P = 5.7-6.35. The error bars for the circles denote the statistical uncertainty of We also plot Klassen’s data at 

mo = mo*'”“ for /3 = 5.5-5.8 as open diamonds. 
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FIG. 10. Results of ro/us. The left figure shows typical fit range (rmin) dependence of ro/fls and its averaged value. The 
right is the result of Os/ro as a function of /3 and its fit curve Eq. (|5E|). 




FIG. 11. The left-hand side shows /3-dependence of the lattice spacing. The solid line is the fit curve Eq. ([H^, while dotted 
and dashed lines are spline interpolations to square and triangle symbols respectively. On the right-hand side /oTJ^ and 

as a function of are plotted. 
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FIG. 12. Comparison of results for various (Miat,C) tunings. The scale is set by ro. The data points are slightly shifted 
along the horizontal axis for the distinguishability. 
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FIG. 14. S-state hyperfine splitting AM(l®5'i — l^So). Results obtained with ci, = UsCs = 1 are collected here. Our results 
are shown by filled symbols for each input, while results by Klassen (set B) and Chen (set C) with the ro input are shown by 
open symbols. In the legend, we give the choice of the anisotropy (( tuning, tadpole factor and scale input. These captions 
also apply to the figures that follow. 



FIG. 15. P-state fine structure splitting AM(l®Pi — l^Po). 
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FIG. 16. P-state fine structure splitting AM(1®P2 — l^Pi). 
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FIG. 19. Spin averaged IP — 1.5 splitting. In the figures, we always omit the bar for the spin average. 
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FIG. 20. Spin dependent IP — IS splittings. 
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FIG. 22. Spin dependent 2S — 15 splittings. 











FIG. 23. Spin dependent 2P — IP splittings. 
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FIG. 24. Charmonium spectrum in the continuum limit. The scale is set by IP — IS splitting. 
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FIG. 25. Rhfs with Cs = Ifu and ^ at ^ = 3 and 2. The thick symbols are the results with = mi, while the thin 

symbols are those with rriq = m 2 . The results with Cs = Ijv but ^ (where m, = mi = m 2 ) are also shown by the dotted 

line = 3) and dashed line (^ = 2). 
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FIG. 26. Klassen’s results of S-state hyperfine splitting AM{l^Si — l^S'o) with Cs = Ijv and C, = (set D). The scale is 
set by vq. Lines denote a^-linear extrapolations. 
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FIG. 28. The results of S-state hyperfine splitting AM{l^Si — 1^5'o) with Cs = 1. The scale is set by ro. 





